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ON HYBRID AND MIXED FINITE ELE:£EKT METHODS 


Th. H. H. lian 

Massachusetts Institute of Technology 
Cambridge, Massachusetts, L'.S.A. 


Three versions of the assumed stress hy- 
brid model in finite element methods and the 

V 

corresponding variational principles for the ' 
formulation are presented. Examples of rank de-^ 
ficiency for stiffness matrices by the hybrid 
stress model are given and their corresponding 
kinematic deformation nodes are identified, A 
discussion of the derivation of general semi- 
Loof elements for plates and shells by the hy- 
brid stress method is given. It is shown that 
the equilibrium model by Praeijs de Veubeke can 
be derived by the approach of the hybrid stress 
model as a special case of senii-Loof elements. 



1, INTRODUCTION 


It is generally realized that an assumed stress 
hybrid element is based on the complernentai-y energy 
principle using equilibrating stress field in Lhe inte- 
rior of the element and independent displacements along 
the element boundary. However, when a compatible dis- 
placement field can be constructed it is also possible 
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to use the Hellinjjer-Heissner principle for the for- 
nulation. Indeed, it may not be Generally knov/n * ^ ^ that 
the original derivation of the assumed stress hybrid 
element v.as made by using the Kelli nger-Reissner prin- 
ciple. In that derivation the assumed stresses happened 
to satisfy the equilibrium equations and the resulting 
element stiffness matrix is identical to that by the 
complementary energy principle, i. e. by a model which 
was later named hybrid model. Hybrid and mixed models, 
thus, are not mutually exclusive. 

In this paper we discuss several problems associa- 
ted with this hybrid/mixed method: 

1) Different formulations of the assumed stress 
hybrid/mixed elements. 

2) Rank deficiency in the assumed stress hybrid 
elements. 

Semi-Loof elements for plates and shells by the 
hybrid formulation and formulation of the equi- 
librium element by the procedure of the hybrid 
stress model. 


2. FGRV.Di.ATir/ES OF TIiS HYBRI U/ MIXED ELEMENTS 


Three different versions of the hybrid element 
formulation are presented here. The first one is by the 
Hellinger-Reissner principle which can be expressed as: 

= StdfioDdry (1) 

where 

0 = stresses 
u = displacements 
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S = elastic compliance matrix 

D = matrix of differential operators that defines 

the strain displacement relations c = D u 

T = vector of surface tractions 

V = volume of the nth element 
n 

= boundary of the nth element over which trac- 
^ tions T are prescribed 

Here $ denotes sxmmation over all elements, and for 
n 

simplicity in the present illustration, body forces are 
considered absent and displacements u a'"e assumed to 
satisfy the prescribed boundary conditions. 

V.hen the state of stress a is in equilibrium, i.e. 


D 0 = 0 in Vp : and fg=T on dVp 


T 

where D represents matrxx of differential operators 
and ^ represents the directional cosine of the surface 
normal, then by the diver£,ence theorem the variational 
functional becomes 

^ J ( jk f b £ I'y dS fu d5]=- (3) 


where dV^ is the entire boundary of the nth element, 

the displacements along the boundary are now devot- 
ed by u. It is noted that the above expression would 

be identical to -r where is the functional cor- 

mc me 

responding to the modified complementary energy prin-^ 
ciple associated with the assumed stress hybrid raodei^J 
When the body force term is included, a mixed varia- 
tional principle results. This also suggests that body 
force can be distributed rationally based on independ- 
ently assumed displacement functions. 



' T'. *' - •' 


/ • 

.till 



I 


In case that both the equilibrium and compat- 
ibility conditions are satisfied both stresses and dis- 
placements can be expressed in terms of the same func- 
tions, The volume integral in Eq. (3) can be reduced 
into a surface integral along the boundary and another 
modified complementary energy principle can be written 


as follows 
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=stntionury (4; 

In the above expression the tractions T and displace- 
ments u arrive from the same stresses a 

The steps to be taken in formulating the assumed 
stress hybrid elements by these three variational prin- 
ciples are indicated in Table 1 . Here the approximation 
of the stresses in the interior of the element is a 
common step for all methods, and although many terms in 
these variational functional are different, the result- 
ing TT -functions are of the same form, i.e. in terms of 
stress parameters fi and nodal ’displacements q. Since 
the stresses are independent for different elements, 
one can take variation with respect to the s in the 
element level and obtain expressions of g in terms of 
q. The final expression for n is now in terms of only 
nodal displacements q as unknowns. The matrix k is the 
element stiffness matrix given by 

k-6^H''G (S) 

^ ^ 

Matrices k can be assembled to form the stiffness ma- 
trix K for the global system. Teking variation of n with 
respect to unrestrained nodal displacements leads to 
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the system of matrix equations of the finite element 
method: 

^3=Q ( 6 ) 


An investigation of the relative efficiency be- 
tween the formulations of the element stiffness matri- 
ces by TTo and for rectangular block elements with 

h me r >1 ^ 


8-nodes was conducted 
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and the results indicated 


clearly that the computing time required for the -n^ 

formulation is shorter than that by the n formula- 

mc 

tion. The reason for this is that in the evaluation of 
each element of the matrix G in Table 1 , the formula- 
tion by involves a sin£;le volume integral while for 
the formulation of the corresponding surface inte- 

gral must be divided into six separate ones for the six 
individual faces. 

The formulation of the assumed stress hybrid ele- 
ment by TT^ has also been extended to the general 6-node 
hexahedron elements which have straight edges but may 
have non-flat faces. In formulating such a finite ele- 
ment by it is required to introduce over the ele- 

ment a set of curvilinear coordinates tj , ^ in the 
same way as the isoparametric element in the conven- 
tional assumed displacement method.. It is v/ell known 
that for conventional isoparametric elements, it is, in 
general, not possible to evaluate, in closed-form, the 
integrals required in the element stiffness matrix, 
hence numerical integration procedures are required. 
However, in the formulation of such elements by the 
assumed stress hybrid models the matrix H can be eval- 
uated in closed form. By recognizing identical elements 
in the H matrix, it is possible to optimize the com- 
puting effort and as a result, the computer time needed 
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to evaluate the .itiffneaa r.iatrix by is only a few 
percent higher than that for the conventional iso- 
parametric element Discussions on the derivation 

of hybrid stress elements by have also been given in 

A typical application of the variational function 

TT ♦ is in two-dimensional linear fracture mechanics, 
me 

Standard technique of complex stress functions in plane 
elasticity problems permits a proper approximate solu- 
tion which satisfies the equilibrium and compatibility 
equations as well as the stress free boundary condi- 
tions at the surface of the crack, A super-element 
which contains an embedded crack can thus be derived 
to be used jointly with conventional finite elements 
for the analysis of elastic stress intensity factors 


3. RAhTC D2FIGIEKCY IR ASSUI/iED STRESS HYBRID ELEMENTS 


In Table 1 the relation between the stress para- 
meters ^and the modal displacement q is governed by 

K i - Gq (7) 

To 121 

It has been pointed out*’ * ^ that if m is the number 

of stress parameters and n is the number of generalized 
displacements of which r nodal displacements must be 
restrained to prevent rigid body motion, then when m is 
smaller than (n-r), kinematic deformation modes v.'ill 
appear and the rank of the stiffness matrix will be 
less than (n-r). It is well known that the condition; 

m > (n-r) (8) 


is only a necessary condition for preventing any kine- 
matic deformation modes, and an element may still be rank 
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Formulation of assumed 
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deficient even if the above inequality is satisfied. 

A kinematic mode corresponds to nodal and bound- 
ary displacements for which no work is done by the as- 
sumed stress distribution, hence is also called zero 
energy mode. This means tho.t 

j5 = h~'g q =0 or Gq = 0 (9) 


Since the assumed stresses are equilibratinc, all ri^id 
body modes will involve no external work. A kinematic 
deformation mode refers to element boundary displace- 
ments v/hich indicate element deformation but involve 
no work from the assumed stresses. 

One of the cases that kinematic deformation modes 
appear when the condition is satisfied is a twelve de- 
cree of freedom rectangular plate element under Kir- 

n 3 1 

chhoff assumption , It is derived by assuminj^ linear 
distributions in stress couples v/ithin 

the element, end cubic distribution in lateral dis- 
placement (w) and linear distribution in normal slopes 

(w ) alone each edge. The condition of Eq. (6) is sat- 
» ^ 

isfied. The element has three ri^id body degrees of 
freedom yet the stiffness matrix has five zero eigen- 
values. One can verify that the boundary displav-.ements 
depicted in Pig. 1 and represented by the following 
equations will yield no work due to linear distribu- 
tions of moments M and M and hence uniform distribu- 

n ns 

tion in Kirchhoff shear along the boundary: 

diong y = b, w = -x(d^-x^) i W.y = 0 
y -~b. w=x(a^-x^); W,y=0 

3long X w = 0; Wy=2a^y/b 

X=-a, w = 0i w^)^-2a^y/b 
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Another kinematic mode Involves similar displacements 
with nonzero values of w eind w along x >± a, and y ■ 
rb respectively. To suppress these kinematic modes it 
is only necessary to add stress terms such as =^^^xy 
and or ''''^ 13 ^^ which v/ill lead 


y ^11*' xy ^12 
to linear distribution of 
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the midpoint of each side as nodal displacements ^ 
aa shown in Pig. 2. For the twelve nodal displacements 
shov/n in the figure, the distribution of normal rota- 
tion w along all edces are in the form of 

» 


% = f(s> = cn-6(s/f) + 6(i/#yj 




=-j-sina, 

C? 

- — cos a, 


w = 0 ^ = -T 

I I 

=C, 


w=0 
^x =^i 
^y 

C, . Cz 
= -ySin a2-y COS dz 


w = 0 

^x = 

W, = c, 


Pi^r, 2 Kinemc-tic defori.xi ticn modes of a 12- 
LO? plate element under linear moment 
distribution 





where s is the coordinate along the edge and t is the 
length of the edge. Here f(s) is symmetric about the 
midside point and its integral over ^ vanishes, thus 
the integral over the product of f(s) and any linear 
function in s v;ill be zero. Also, along all edges the w 
distribution is cubic ond is antisyi.jnetric about their 




midpoints. Now under linear distributions in stress 
couples, the bending moment is linear along s and 
the Kirchhoff shear is constant. They should do no 
work under the boundary of displacements described 
above. Since w is zero at all corners, the corner 
forces also do no work. The deformation pattern shown 
in i^ig. 2 represents the combination of two independent 
kinematic modes and the rank of the stiffness matrix is 
seven. 

Finally, a remark should be made about the rectan- 
gular membrane element derived by using 5 P terms which 
was used as an example when the assumed stress hybrid 
element was first introduced. The element has eight 
degrees of freedom and three rigid body modes, hence, 
Eq. (bj IS satisfied and for general arrangement of the 
reference axes the resulting element stiffness matrix 
will have a rank of five, iiov/ever, for a square element 
if the diagonal lines are used as the reference axes as 
shown in Fig, 3, the resulting stiffness matrix will 
have a rank of only three. If the stress assumption is 

=^, ^ $4y 

Oy^Bi * (12) 

r^y^Bj 

then in the - n coordinate system the corresponding 
assumed sti'ess is 

Orr 02*04 as n (13) 

Tjfy = Qj -Q4rj-as i 
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i'ig. 3 Squui'O n:e.7.brtine element with reference 
txes for stresses tilong the diu^-onals 

In this case tne following two nodal displacement pat- 
terijc are kinematic modes: 


Ca) 
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- -% = % = 

=1 


and 
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= % = = 

% =0 

(b) 
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= % = % = 
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Similarly, if the x-cxis coincidec with one diagonal of ' 

a general rectangular membrane element, the resulting * 

element stiffness matrix by 5-^ assumption will hove | 

Q rank of only four. * 

Spilker et have studied the problem cf 

rank deficiency ci the 5-^ square element by first 
identifying all independent deformation modes that can 
uniquely define any displacements of the element and 
then substituting into Eq. (9) to see whether G, q 
vanishes. By changing the angle 0 between x and C axes 


they werj able to detect the tv/o kinematic deformation 
at 6 > 43 « This example cn t. e membrane element also 
Illustrates the danger of not takinij complete poly- 
nomials in the stress terms in the hybrid element for- 
mulation. In that case, the resultinc element stiff- 
ness matrix v/ill be lack of invariance on and, indeed, 
it may be rank deficiency under certain reference coor- 
dinates. 

For any newly developed hybrid stress elements it 
is succ®sted, in addition to latch- test, an eijenvalu* 
survey should be conducted in order to detect any ki- 
nematic deformation modes. If invariance of the stil?'- 
noss matrix cani;ot Le mainta .led, such a survey will 
have to be made for different coordinate system.s used 
for the stresses. It should be noted that kinematic 
deformat: in modes can always be suppressed by adding 
appropriate stress terms. 

4. £a4l-LC0F F: EI>i:i;T3 FOR PLATES ALT) SHELLS 
BY THE HYBRID F0fU4ULATI0N 

For the analyses of plates and shells, one of the 
difficult tasks is to match the compatibility at a node 
at which the reference planes of the elements are not 
coplanar. At such nodes ell six degrees of freedom 
should be considered but for plate and shell elements 
a maximum of five degrees of freedom can be used at a 
node. Irons^^^^ has suggested the use of the so-called 
semJ.-Loof element for which the normal rotations along 
each edge are defined at nodes which are not located at 
corners of the element. By the conventional assumed 
displacement method it is still a difficult task to 
construct shape functions for a semi-Loof element. But, 




it is a more or less routine procedure to formula’’, e \ •: 

’ I 

such an element for plates and for simple shells by the , 

I 

assumed stress hybrid method. 

The equilibrium model of Praeijs de Veubeke^^^^ I 

always results to semi-Loof elements. It has been shown 
that this model and the assumed stress hybrid model can 
both be derived by the modified complementary energy 

r 2 ) 

principle . The only difference is that for the 
former, tractions along each boundary are represented 
by generalized loading parameters and the corresponding 
fictitious nodal displacements are v/eighted integrals 
of the boundary displacements, while for the latter, 
boundary displacements are interpolated in terms of 
nodal displacements and the corresponding nodal forces 
ere obtained from the variational sense. The expression 
for stiffness matrices are the some for both models ex- 
cept the G matrices in £q. (5) are derived differently 
for the t;vo models. 

.".e can now show that the equilibrium model can 
also be formulated by boundary displacement interpola- 
tion in the same way as for tt in Table 1 . Consider 

me 

a problem for which certain nodal displacements are of 
the Semi-Loop type, hence, the corresponding boundary 
displacements are independent from one edge to the 
other. This means that for each element the boundary 
displacement continuity is not maintained at the cor- 
ners. 

For simplicity let us consider a straight boundary 
0<s<C, with its traction T (s) represented by a poly- 
nomial of order m. 

Such traction distribution can be represented by 
m+1 nodal forces arbitrary locations 

s^,..., when the following conditions are satisfi- 
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In 


I Qi 5--J'T(s)s''ds. forn=0.-, m (I 

1=1 

We can easily show that if the boundary displacement 
u(s) is also a polynomial of order m then 

T(s) u(s) ds = T Qj Uj (I 


where is the nodal displacement at s^. Conversely, 
if the boundary displacement is interpolated timroutjh 
the nodal displacements at the ra-^1 nodes alone the 
boundary, the corresponding nodal forces obtained by 
Sq. (15) will satisfy the conditions ^iven by iiq. (14), 
hence, v/ill be equipollent to the boundary traction 
distribution T(s). Ahen the boundaries of two nei~h- 
borin^ elements are approximated by polynomials of the 
same order, a compatibility of the nodal displacements 
will then necessarily guarantee the reciprocity of the 
nodal forces and, hence, the pointwise interelement 
equilibrium condition is maintained. The equilibrium 
model by Praeijs de Veubeke, thus, is a special semi- 
Loof element by the assumed stress hybrid model when . 
the number of nodes along a boundary exceeds by one the 
order of the polynondnal which represents the corre-* 
spending tract' on. An exception in thi" case is that 
for constant traction distribution there is only one 
node v/hich should be located at the midpoint of the 
edge. For example, for the linear moment triangular 
plate element by Praeijs de Veubeke^^*^^ , the nodal dis- 
placements are lateral displacements w at the three 
corners and the three midside points and two normal 







rotations alone each ed^e. Here, the normal bending 
moment is linear along each edge, hence, by using 
tv;o corresponding nodal displace, ents w _ along the 

9 ^ 

edge, the equilibrium conditions are satisfied exactly. 
It should be remarked that in this case the locations 
of the semi-Loof nodes may be arbitrary. 

Por the equilibrium model, kinematic deformation 
modes often exist in one element or in a group of ele- 
ments and appropriate superelement technique is often 

[ ^ 0 ] 

required to suppi-ess the mechanisms . On the other 
hand, for the same choice of nodal displacements, the 
number of stress terms can be increased by the assumed 
stress hybrid model, hence, the possibility of any 
kinematic deformation modes can always be eliminated. 
Thus, the assumed stress hybrid method is a most prom- 
ising approach for the seni-Loof element. 


5. COHCLLTIKG RiyjuulS 


Por the formulation of the assumed stress hybrid 
finite elements there are different methods of approach 
to choose from end the selection of the appropriate 
stress fields requires sufficient physical insight. It 
is not likely that, in the future, hybrid stress ele- 
ments can be formulated as simply as the use shape 
function routines and numerical integration techniques 
in the conventional assumed displacement method. But 
abundant knov/ledge about this method has been collected 
and future development of hybrid stress elements should 
be and can be made fool proof and should be implement- 
ed according to the moat efficient method of approach. 
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